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Gravitational Microlens Mapping of a Quasar Accretion Disk 
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Abstract 

Since the radiation from different portions in the central region of a quasar can be successively amplified 
during a microlensing event, microlensing light curves provide fruitful information regarding the emissivity 
distribution of an accretion disk located at the quasar center. We present a basic methodology of how 
to map the emissivity distribution of the disk as a function of the radial distance from the center, Q(r), 
from 'observed' microlens light curves during a caustic crossing event. Our proposed method is based on 
the standard inversion technique, the so-called regularization method, and Abel's transformation of the 
one-dimensional luminosity profile integrated along the line parallel to the caustics. This technique will be 
used to map the disk structure in Q2237+0305, for which the HST and AXAF observations are scheduled. 
A reconstruction of the image on length scales of several to ten AUs is quite feasible for this source, as long 
as the measuring errors are within 0.02 mag and the observation time intervals are one week or less. 

Key words: Accretion disks — Active Galactic Nuclei — Black Holes — Microlensing — Quasars: 
individuals (Q2237+0305) 



1. Introduction 

In many cases of astrophysical objects, their direct 
images are too small to resolve with the usual tele- 
scopes. In the case of a quasar accretion disk, for ex- 
ample, its angular size is, typically, on the order of 
#d ~ O.Olpc/lGpc = 10~ n rad « 1/xas. In the case of 
a binary accretion disk, on the other hand, we estimate 
8 d ~ 3 x 10 10 cm/lkpc= 10~ n rad« 1/Ltas. These are both 
considerably below the angular resolution of any present- 
day telescopes. What is usually attempted is, therefore, 
to construct theoretical models based on the basic equa- 
tions of (magneto)hydronamics and to infer their struc- 
ture from the limited number of information, such as the 
total radiation output originating from their entire sur- 
face at some wavelength bands. 

In some special cases, however, we can 'resolve' the 
spatial structure of the disk by using a 'natural' tele- 
scope. A good example is the technique of eclipse map- 
ping (Home 1985). Since a certain fraction of an accre- 
tion disk is shadowed by a companion star in an eclips- 
ing close binary and since the shadowed part varies with 
time in accordance with the orbital phase, we can map 
the disk emissivity distribution from the eclipsing light 
curves. A variety of disk luminosity profiles has been re- 
vealed with this technique and our knowledge about the 
disk structure has been remarkably enriched. Unfortu- 
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nately, however, this technique is irrelevant to probing 
a quasar accretion-disk structure because there are no 
eclipsing objects known to date (see, however, McKer- 
nan, Yaqoob 1998). 

An alternative and potentially powerful method using 
a microlensing phenomenon has been considered. Chang 
and Refsdal (1979, 1984) considered flux changes that 
might occur due to gravitational lcnsing by a single star 
in an extended gravitational lens galaxy. Now, a stim- 
ulating discussion has started that such a microlensing 
phenomenon can be used to investigate the quasar central 
structure (Grieger et al. 1988; Rauch, Blandford 1991) as 
a 'gravitational telescope' (Blandford, Hogg 1995). The 
Einstein Cross, Q2237+0305 (e.g., Huchraet al. 1985), is 
the first object in which quasar microlensing events were 
detected (Corrigan et al. 1991, Houde, Racine 1994, see 
also Ostensen et al. 1996). These observations suggest 
that microlensing events seem to take place almost every 
year. This rather high frequency is consistent with the 
microlens optical depth of r ~ 0.2 — 0.8 obtained by a re- 
alistic simulation of the lensing galaxy (i.e., Wambsganss, 
Paczyhski 1994). 

We, here, specifically consider the microlensing event 
of this source caused by the so-called 'caustic crossings.' 
Several authors have already calculated and inspired this 
'caustic' case based on simple disk models (e.g., Wamb- 
sganss, Paczyhski 1991; Jaroszyhski et al. 1992) or re- 
cently on the realistic disk models (Yonehara et al. 1998, 
1999). All of these calculations correspond to the so- 
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called forward problem; i.e., they calculated microlens 
light curves based on disk models prescribed a priori. In 
contrast, we are, in the present study, concerned with 
a distinct approach called the inverse problem; i.e., we 
consider how to reconstruct the disk image (or more pre- 
cisely, the emissivity distribution) of a quasar accretion 
disk from microlcnsing light curves. 

For this purpose, we first simply apply the regulariza- 
tion method, one of the most well-known non-classical 
inversion techniques. The basic methodology of this 
method was already given by Grieger et al. (1991, here- 
after GKS). Indeed, their proposed procedure has been 
very successful, but they have assumed rather idealized 
situations and, thus, there remain some problems requir- 
ing further investigation before realistic applications are 
obtained. 

1. GKS assumed no radiation from any part outside a 
caustic. As a result, they could precisely determine 
the contribution from a small outer portion of the 
disk to the total light from the shape of the light 
curve at the beginning of a microlens event. If this 
were the case, we would not be able to see the disk 
emission in the absence of a microlcnsing, contrary 
to the observations. 

2. The goal of GKS was to reconstruct a one- 
dimensional luminosity profile [P(£) in their nota- 
tion] which is the luminosity integrated over the disk 
plane in the direction parallel to that of a caustic (see 
figure 1 for the definition of £). To make a direct 
comparison with accretion disk models, however, it 
will be more convenient to express the emissivity 
distribution as a function of the distance from the 
center, Q{r), since the central part of the disk can be 
reasonably assumed to be axisymmctric. It is thus 
needed to transform P(£) to Q(r). 

3. GKS assumed a rather smooth emissivity distribu- 
tion around the center, which certainly makes the 
analysis easier than otherwise. According to realis- 
tic disk models, however, it seems more likely that 
the emissivity has a power-law dependence on the ra- 
dius, thus being sharply peaked around the center. 
We wish to know what fraction of the total energy 
output originates from a compact region in the disk. 

The aim of the present study was to improve the GKS 
method so that it can be applied to more realistic sit- 
uations. Our version of the inverse technique will be 
described in section 2. We then present the results of 
disk mapping for the calculated microlcnsing light curve 
in section 3. The final section will be devoted to a sum- 
mary and discussion. 

2. Inversion Procedures 
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Fig. 1.. Schematic view of a 'caustic' crossing during a mi- 
crolens event. The caustic is represented by the thick vertical 
line, and the left-hand region corresponds to the parts inside 
the caustic and is subject to microlcnsing light amplification. 
One-dimensional luminosity profiles, Pi, P%, ■ ■ ■ , Pj, are the in- 
tegral of the local emissivity of the disk along the vertical lines 
parallel to the caustic. 



2. 1 . Setting Matrix Forms 

We first describe the regularization method adopted 
by GKS with some modifications. 

Chang and Refsdal (1979, 1984) have shown that the 
flux from the bright double image of a point source close 
to a caustic can be amplified approximately according to 
(x — £) _ j wnere x an d £ are, respectively, the radial 
distances of the caustic and that of a part of the disk in 
question from a line which is parallel to the caustic and 
which crosses the disk center (see figure 1). Note that x 
is time-dependent; i.e., 



x = x(t) = V caus £/r , 



(1) 



where the unit of length, ro, is taken to be a typical disk 
dimension, t = corresponds to the time when a caustic 
crosses the center of the disk, and V^ aus is the transverse 
velocity of the caustic, also including that of the peculiar 
motion of the foreground galaxy relative to the source 
and the observer, on the disk plane. The total observed 
flux will then be approximated by 



F(t) = F f 



A( X - omn. 



(2) 



Here, Fq is a constant representing the total disk flux 



outside the caustic (x < £); i.e. 
amplification factor is given by 

A(x - = ' 



JP(0d£, = 1, and the 
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where k denotes the amplification factor inside the caus- 
tic (a; > £) and depends on the spatial distribution of 
lcnsing stars and the number distribution of lens masses. 
Note that GKS set A = outside the caustics and, hence, 
k is arbitrary, which grossly simplified the analysis. Also 
note that inclination effects arc included in the expression 
of F and V caus (discussed later) . 

We assume that P(£) is represented by a continuous, 
piecewise linear function, 



P(0 



& - 0-i 



P.i-1 



£-&•- 



t±p j 

- 0-i 



(4) 



for 6_i < £ < £, 



For simplicity, we take equal-distant mesh points normal- 
ized by the disk radius, 



C? = 1,2,---,J), 



(5) 



with A£ = 2/ J, and we take the total number of mesh 
points (J) to be equal to that of the observed time se- 
quences (£,;, i = 1, 2, • • • , J). If we set Pj = F(ti), equa- 
tion pi) leads 



F = F Q KP or F i = F }^K ij P 

3=1 

(i=l,2,.--,J) 



(6) 



with 

P=(P 1 ,F 2 ,-.., J F J ) T , P=(F 1 ,P 2 ,-..,P J ) T , (7) 
and K being a J x J matrix given in Appendix. 

2.2. Regularization Technique 

It seems straightforward to derive the emissivity dis- 
tribution, P(£) or Pj (J = 1, 2, • • • , J), from equation (ra) 
by calculating the inverse matrix of K, 

'-HF/F ). (8) 



P = K" 



However, this does not work efficiently, since the observa- 
tional data usually contain measuring errors, which will 
be greatly amplified when calculating the inverse matrix 
(GKS). 

To resolve this issue, a regularization technique was 
proposed. First, we express the observed flux as 

j 
F = P KP + 5 or F t = F ^ KyPj + S t , (9) 

i=i 

with Si being the error in measuring the flux, Pj. Next, 
as a measure to evaluate how the P(£) profile is smooth 
we introduce a badness function, 



L(P) = J 



d 2 P(0 



d£, 2 



r/e 



(10) 



(Note that L(P) = if P(C) is a linear function of £.) 

Then, the problem is reduced to determining the func- 
tional form of P(£) which gives a minimum value of 



[(P/P )-KP] 2 + AL(P). 



(11) 



Namely, there are two important factors to be considered: 
fitting to the given light curves (P, the first term) and 
smoothing the emissivity profile (P, the second term), 
where L(P) is the matrix form of L(P), 



L(P) - E 



P 



3+1 



2P + P,-_ ! 



(AO 2 



Ae 



(12) 



The controlling parameter, A, is called the smoothing 
parameter. In the limit of vanishing A, the solution, P(£), 
is given by equation (||) which can precisely reproduce 
the observed light curves but is not always very smooth, 
especially in cases with large measuring errors (GKS). In 
the limit of very large A, on the other hand, the solution 
is smoothest, since it gives a straight line on the [£, P(£)] 
plane, but may not give an excellent fit to the observed 
flux variations. We thus need to choose a moderate value 
of A which satisfies 



1 J 

2—1 




J 



3=1 



Fo^KijPj 



(13) 



In other words, the final solutions should be as smooth 
as possible under the constraint that a reproduced flux 
variation agrees with the observed variation within the 
error bars. 

After some algebra, minimizing equation (llll) is equiv- 
alent to setting 



(K T K + AH)P = K T (P/P ), 



(14) 



where H is a matrix (see Appendix of GKS). The final so- 
lution, P(£), is obtained successively by solving equation 
jl4| ) for a given A satisfying equation (|l3|). 

These procedures are performed for a specific value of 
k, say, k = 0.5 (which gives a variation with an amplitude 
of ~ 0.5 mag for Q(r) ex 1/r). We repeat the same 
procedures with different values of k and determine the 
k- value so that L(P) reaches the minimum, yielding the 
smoothest P(£). 

2.3. Transformation from P(£) to Q(r) 

As stated in Introduction, our goal is to determine the 
emissivity distribution of the disk as a function of r (not 
£), where r is the radial distance from the center of the 
accretion disk in the unit of tq (r = 1 at the outer rim) . 
Thus, we need to transform P(£) to Q(r). 
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Note that the dimension of Q(r) is flux; i.e., the energy 
emitted per unit time from a unit surface area. Then, 
P(£) can be expressed as an integral of Q(r) as 



^(0 



Q(r)drj 



Q(r) 



rdr 



v^e 



(15) 



It is of great importance to note that equation (llE 
takes the form of Abel's integral; that is, the in- 
verse transformation is straightforward (see, e.g., Binney, 
Tremaine 1987). Since we have two independent sets of 
P(£), one at £ < and the other at £ > 0, we can sepa- 
rately obtain two sets of Q(r), Q~(r) and Q + (r); namely, 



Q~(r) 



'(r) 



dP(Q d£ 

dP(0 



<% 



and 
d$ 



ve 



(16) 



Note that quite generally dP/d£ > for £ < and 
dP/d£, < for £ > 0. 

The next procedure is used to derive an expression for 
)±(r fe ) in terms of P, = Pfo) (j = 1, 2, • • • , J) 



>t 



with rfc = (2k 
Q~(r k ) = 



l)/Jfor k= 1,2, 

(J/2)-* 



• , J/2. We have 



IT £-~t 



P 



J + l 



£ 



(1/2) 



1 
+ - 



Pi 



and 



Q + (r k ) 



Tryr - 



J-i 

£ 

j=(J/2)+k 
1 P/ 



(17) 



P, 



■i+i 



V^ 2 +(l/2) 



tt^/t 



(18) 



with £j+(i/ 2 
used £i/ 2 = 



find Q (rfe)'s for k 
(j = l,2,---,J). 



= Ife + &+i)/2| = (2j - J)/ J, and we 
-1 and Cj+(i/2) = +1- We can successively 



l,2,---,(J/2) - 1 from Pfo)'s 



£.^. Producing Light Curves 

To test how the above procedures work, we calculate 
the expected flux variations based on specific models for 
Q(r) [hereafter denoted as QmodeK^) to distinguish from 
the reconstructed values, Q (r)]. It is important to note 
that accretion disk structure can well be described as be- 
ing self-similar (or more precisely, each physical quantity 
is expressed as a power-law function of the radius, see, 
e.g., Shakura, Sunyaev, 1973; Narayan, Yi 1995). We, 
hence, prescribe the emissivity distribution as 



cdelO") oc r~ 



for 0.01 < r < 1.0 



(19) 



with a being a positive constant. Note that a — 3 for 
the standard- type disks, whereas QmodeiM is much flat- 
ter, aSil, in optically-thin, advection-dominated disks 
(Manmoto et al. 1997). 

Light curves are calculated by 

F(U) = Po / rdr d0A(x - r cos0)Q modcl (r) 
Jo Jo 

+S % (20) 

with the amplification factor [A(x — rcosO)] being given 
by equation (]3j). We continuously change the angular 
separation between the caustic and the center of the ac- 
cretion disk (see Eq.[l]]). The unit of time (to) is taken 
to be the crossing time over which the caustic moves 
on the quasar disk plane from £ = to +1. That is, 
to = ^o/K;aus, and, hence, x(t) = t/to from equation 
(H). On the intrinsic light curve, we add measuring er- 
rors, Si. We calculate three models for each QmodeK?*) 
prescription; F(t) at each time is randomly fluctuated 
around the mean value in a Gaussian way with assumed 
standard deviations, Am (in magnitude). 

We plot in figure 2 how different microlens light curves 
are produced by changes of the QmodeK^) prescription. 
The normalizations are taken so as to give P = 1.5Pq at 
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Fig. 2.. Microlens light curves of disks whose emissivity profiles 
are Q(r) oc r~ a with a = 0, 1, 2, and 3. The values of a are 
indicated in the figure. The larger is a, the more sharply peaked 
is the light curve. Note that the time of the peak flux is t = 
when a > f , but it is shifted to t ~ 0.7*o when a = 0. 

t = O.Uq. Importantly, the peak luminosities are reached 
at the time of caustic crossing over the disk center (t = 0) 
only for a > 1, while a disk with a flat emissivity profile 
(with a — 0) yields the peak at a later time, t ~ 0.7to- 
This is because for the flat disk the amplified area (which 
reaches its maximum at t = to) is also an important 
factor. For a > 1, in contrast, the emergent flux is rather 
insensitive to the extent of the amplified area, since a 
large fraction of radiation comes from the very center of 
the disk. 
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Am=0.10mag 

i i i i I i i i i I i i i i I i i i i I i i. 



Am=0.05mag 

i i i i I i i i i I i i i i I i i i i I i i. 



Am=0.02mag 

i i i i I i i i i I i i i i I i i i i I i i. 
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Fig. 3.. Results of the reconstruction for the cases with Qmodol °c 1/r. The dashed curves in the lower three panels are the prescribed 
cmissivity distribution (Q mo dcl tx i/ r )- The solid curves in the lower panels at r < (or r > 0) represent the reconstructed Q~ (Q + ) 
profiles calculated from P(£) at £ < (£ > 0) for three cases with Am = 0.1, 0.05, and 0.02 from the left to the right, respectively. 
The middle panels are the prescribed P(£) calculated from QmodelM, an d the reconstructed ones from the light curves above. The 
upper three panels show the original light curves calculated by the prescribed Q mo del( r ) (by the dashed lines) with error bars, and 
the reconstructed ones (by the solid lines). 



To summarize, those cases with a — and a > 1 give 
distinct flux variations. It is thus important to distin- 
guish these two critical cases in disk mapping. 

3. Results of Image Reconstruction 

The results of the reconstruction are displayed in figure 
3 for the case with Q m odc\(r) oc 1/r. The prescribed 
emissivity profile is illustrated by the dashed lines in the 
lower three panels. To calculate the flux variation, we 
set a constant time interval, At = 0.2to, and the total 



number of observations is J = 24. The calculated (i.e., 
'observed') flux at each time is displayed together with 
error bars in the upper three panels by the dashed lines 
for different magnitudes of the mean errors, Am = 0.1 
(left), 0.05 (middle), and 0.02 (right), respectively. The 
dashed lines in the middle ones show the one-dimensional 
emissivity, -P(£), calculated from Qmodei(0- 

Similarly, the reconstructed F(t) and P(£) are shown 
in the upper and middle panels with the solid lines. The 
two solid lines in the lower panels display Q~(r) at r < 
and Q + {r) at r > 0, respectively. Note that the recon- 
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structed fc-value does not always reproduce the original 
value of fc = 0.50; we find k = 0.61, 0.51, and 0.47 for 
Am = 0.10, 0.05, and 0.02, respectively. The recon- 
structed flux variations and emissivity profiles are not 
very sensitive to these small deviations in the fc-values, 
however. 

For a = 1, large errors (Am = 0.1 mag) tend to pro- 
duce a rather flat P(£) profile and thus a smooth Q(r). 
This is because a flat model [P(£) = const.] is compat- 
ible with the light curves within the error bars so that 
the technique prefers a flatter P(£) profile (see figure 2). 
Note again that the peak shifts from t — in the original 
light curve to t ~ 0.7io in the reconstructed one. Such 
a problem does not arise for those cases with smaller er- 
rors. To reproduce a steep Q(r) profile up to the inner 
parts, Am = 0.02 mag is necessary. It depends on the 
observing intervals how close to the origin the mapping 
technique can reproduce the original image. Surely, fre- 
quent observations, especially at around the peak flux, 
are preferable (GKS). 

Figures 4 and 5 display the results of those cases with 
Qmodoi oc 1/r 2 , and ex 1/r 3 , respectively. In each figure, 
we omit the cases with Am = 0.05 and P(£) plots to 
avoid any complications. 
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Fig. 4.. Results of the reconstruction and the original for the cases 
with Q mo dcl (X- l/»* 2 . The lower two panels are the prescribed 
emissivity distribution (by the dashed lines) and the recon- 
structed one (by the solid lines) for two cases with Am = 0.1 
(left) and 0.02 (right), respectively,. The upper two panels show 
the original light curves (by the dashed lines) with error bars, 
and the reconstructed ones (by the solid lines). 

The reconstructed fc-values are 0.63 (Am = 0.10) and 
0.74 (Am = 0.02) for a = 2, and 0.76 (Am = 0.10) 
and 0.78 (Am = 0.02) for a = 3, respectively. These 
all deviate greatly from the original value of k = 0.5, 
particularly when a is large or Am is small. Nevertheless, 
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Fig. 5.. Same as figure 4, but for cases with Q mo dcl °c 1/r 3 . 

the inner steep rise parts are reasonably well reproduced 
by the mapping. The wing parts are poorly reproduced 
for a steep Q(r) profile. In fact, Q(r) sometimes goes 
below zero. This tendency is rather enhanced for small 
Am's. 

To conclude, in cases with sharply peaked emissivity 
profiles we can still obtain reliable information regarding 
the extent of the substantially emitting region, although 
we cannot trust the results about the outer zones sur- 
rounding the central bright zone. 

4. Summary and Discussion 

Let us consider specifically the case of Einstein Cross 
and thus insert the model parameters relevant to this 
source. The Einstein-ring radius on the source plane is 



r E = 0eD o 



2 ) 






A>1 



1/2 



i.sxW^y 72 . 

\M J 



(21) 



where Mi cns is the typical mass of a lens star, and D\ s , 
.Dos, and D \ represent the angular diameter distances 
from lens to source, from observer to source, and from 
observer to lens, respectively. To evaluate these dis- 
tances, we assume the redshifts corresponding to the dis- 
tances from the observer to the quasar and from the 
observer to the lens of, z os = 1.675 and z a \ = 0.039, 
respectively (see Irwin et al. 1989), and also assumed 
an Einstein-de Sitter universe and Hubblc's constant to 
be Hq <~ 60km s _1 Mpc _1 , according to Kundic et al. 
(1997). 
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Another important length is r cross , the caustic crossing 
length over the quasar image plane per observational time 
interval, At; 



v t At— 

2.0 x 10 1; 



v t 



300 km s _1 



n'''' 11 !22i 



where ^t is the transverse velocity of the lens on the lens 
plane (14 aus = v t D 0S /D \). Surprisingly, this is compara- 
ble to the Schwarzschild radius, r g ~ 3 x 10 13 (Af 8 ) _1 cm 
(~ 2AU) for a 10 8 Afg black hole and is much smaller than 
te- Thus, by weekly observations can one determine the 
disk emissivity distributions on length scales of ~ 10AU 
or ~ 5r g (M8) _1 for vt ~ 300 km s _1 . 

To summarize, we have improved the reconstruction 
technique previously developed by GKS in such a way 
that a direct comparison with accretion disk models is 
possible. We have found that for deriving the emissivity 
distribution as a function of r on scales down to several 
to ten AUs, we need an accuracy of Am < 0.02mag and a 
sampling interval within one week. In cases in which the 
emissivity profile is rather centrally peaked (i.e., if a > 1), 
we can still reproduce the inner bright zone reasonably 
well, but cannot trust the results of the outer zone. 

It might be noted that the problem treated here seems 
to be closely related to that of reconstructing stellar 
brightness profile (e.g., limb darkening) from photometry 
of Galactic microlensing events (see Gaudi, Gould, 1999, 
and references therein). The present technique may be 
used to analyze such data. 

The inclination of the disk has two important effects. 
If the disk plane is tilted by an angle of i with respect to 
the line of sight, the apparent disk flux will be reduced 
by a factor of (cosi) -1 . This affects the normalization 
constant, F . Further, if the disk is tilted by an angle of 
iii with respect to the direction of the motion of a caustic, 
the apparent transverse velocity of the caustic will be in- 
creased by a factor of (cos«||) _1 . This affects the caustic 
crossing, V cauB , and thus the time and length scales, to 
and tq. In other words, the normalization constants both 
in the ordinate and abscissa in the (t, F) diagram and 
thus in the (r, Q) diagram are subject to the inclination 
angles. There remain ambiguities in the normalizations 
of V{) and Qq. Since it is difficult to evaluate these in- 
clination angles and other uncertain factors, we rather 
focused our effort on the shape of the non-dimensional 
emissivity distribution, Q/Qo, as a function of r/rg. 

We must consider the fact that real photometric light 
curves will be unevenly sampled due to gaps in the 
observing schedule caused by various factors, including 
weather. Therefore, an even more realistic case could 
be made by making a model of such gaps. To see such 
an effect, we made a light curve in which a certain frac- 
tion, ~ 25%, of the regularly spaced light curve samples 
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Fig. 6.. Same as figure 4, but for the unevenly sampled case. The 
given profile is Q mo dei oc 1/r. 



is removed stochastically. The results for the cases with 
Qmodei of l/ r ar e displayed in figure 6. The reconstructed 
fc-values are 0.48 (Am = 0.10) and 0.415 (Am = 0.02), 
respectively. Because there are no significant changes, as 
long as frequent observations are made around the peak, 
a lack of data (uneven sampling) should not cause any se- 
rious problem. Obviously, however, poor sampling rates 
lead to a poor spatial resolution in the disk mapping (see 
equation (2^). The results would sensitively depend on 
whether or not there are enough observation runs at the 
times around the peak. 

HST and AXAF observations of the Einstein Cross are 
scheduled. If observed with multi- wavelength bands, a 
microlensing event should first clearly resolve the multi- 
wavelength radiation properties of a disk in a distant 
quasar on length scales down to several AUs. 
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Appendix 

Equation (0) can be rewritten as 
J 






(23) 
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where fy is non-zero only for Xi > £j—i and is 
fij = I A(x-S)P(Z)d£ 



for (i,j) = (1,1), (1,2),..- 
to derive expression for K, 
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(J, J), it is straightforward 






Kij = E hl+1 + D 
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(30) 



k 



' (j W P '-> 






by inserting equations (26) and 



A£ 



^ «, (24) 



with A£ = £j —£j_i. We introduce matrices Ey and Dy 
such that 



We then have (for Xi > Cj-i) 



riain(£j,Xi) 



(25) 




W-bfo-Xi)] 



4fc ( 3 ; i -g J -- 1 ) 3 / 2 -[p( a ; i -g J -) ] 3 / 2 
3 A^ 

with s = ^j — £, A'£ = p(^j — Xi) and 

p(3i-&)= ( n i_ ^ ^ ^if* 
PV sj; |_ for ^ < £,-. 

Similarly, we find (for x^ > £j-i) 



(27) 
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2A£ A£ 

4fe(^-^-i) 3/2 -[p(^-0)] 3/2 
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(28) 



with u = £ — £j-i and A"£ = min(A£,iz;i — £j-i)- For 
Xi < Cj'-i both are zero, Sy = Dij = 0. Since we have 
from equation (]25|) 



fij — {Eij+i + DijjPj = KijPj, 



(29) 



